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Aims. We aim to provide an estimate of the minimum temperature of the quiet solar chromosphere 

Methods. We perform a 2D radiation-MHD simulation spanning the upper convection zone to the lower corona. The simulation in- 
cludes non-LTE radiative transfer and an equation-of-state that includes non-equilibrium ionization of hydrogen and non-equilibrium 
H 2 molecule formation. We analyze the reliability of the various assumptions made in our model in order to assess the realism of the 
simulation. 

Results. Our simulation contains pockets of cool gas with down to 1660 K from 1 Mm up to 3.2 Mm height. It overestimates the 
radiative heating, and contains non-physical heating below 1660 K. Therefore we conclude that cool pockets in the quiet solar chro- 
mosphere might have even lower temperatures than in the simulation, provided that there exist areas in the chromosphere without 
significant magnetic heating. We suggest off-limb molecular spectroscopy to look for such cool pockets and 3D simulations including 
a local dynamo and a magnetic carpet to investigate Joule heating in the quiet chromosphere. 
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1. Introduction 

The internetwork solar chromosphere is continuously pervaded 
by acoustic waves and shocks with periods of around 3 min- 
utes. A large literatu r e exis ts on this subject. Some exam- 
ples are: iKriiger et al.l (l200ll) . who studied upper-photospheric 
and low-chromospheric oscillations observed in the UV contin- 
uum around 150 nm with the Tr a nsition Region and Coronal 
Explorer (TRACE). iRutten et all (120081) studied the response 
of chromospheric diagnostics to photospheric events such as 
granular buffeting and exploding granules as observ ed with the 
Dutch Open Telescope (DOT). Somewhat longer ago. lLites et alJ 
( 1 1993b studied solar chromospheric oscillations using spectro- 
grams in Call H obtained with the Vacuum T ower Telescope 
at NSO/Sacramento Peak. IWikst0l et all d2000h studied upper- 
chromospheric and transition region oscillations using spectral 
time series obtained around 103 nm with the SUMER spectro- 
graph aboard SOHO. These and other studies confirm the picture 
of acoustic wave excitation by granular dynamics in the photo- 
sphere, the waves then propagate predominantly in the vertical 
direction and evolve into shocks. Depending on the magnetic 
field topology they occasion ally retain their identity all the way 
up into the transition zone dMcIntosh et al.l 1200 ll) . The shock 
fronts are high temperature disturbances in an otherwise cool 
background atmosphere. 

Observatio ns of CO lines in quie t sun areas (e.g., 

Noves & HalJ Il972t lAvres & Testermanl Il981t lAvres et all 
20061) indicate temperatures low enough to form significant 
amounts of CO molecules. Static ID semi-empirical models de- 
signed to reprod uce such observation s can have temperatures 
down to 2750 K (lAvres & RabinllT9 96). These low temperatures 
in such static models correspond to the cool inter-shock phases 
of the quiet chromosphere. 



This scena rio was confirmed in th e one-dimensional (ID) 
simulations of ICarlsson & Steinl (1 19971) that explained the for- 
mation of Can K2V and H2V bright grains as an effect of 
upward-propagating shock waves excited in the photosphere. 
The first attempt to model such wave-response in the chro- 
mospher e with multidimension al numerical simulations was 
made by ISkartlien et al.l (l2000l) . They studied the response of 
the chromosphere to collapsing granules using a 3D radiation- 
hydrodynamic (RHD) simulation spanning from the upper con- 
vection zone to 1.2 Mm above (T500) = 1, and again confirmed 
the picture of the internetwork chromospher e as a cool back- 
ground state pervaded by waves and shocks. IWedemever et al.l 
(120041) performed a 3D RHD simulation with higher spatial res- 
olution but less sophisticated radiative transfer. The higher res- 
olution allowed them to distinguish different shock front shapes 
depending on the excitation mechanism: spherical fronts excited 
by pressure fluctuations in intergranular lanes; planar fronts, ex- 
cited by simulation-box oscillations, the simulated counterpart 
of solar /7-modes; and the previously described irregular front 
shapes caused by collapsing granules. The simulations employed 
simplified radiative transfer and an equation of state (EOS) based 
on Saha ionization equilibrium for hydr ogen, an assumption tha t 
fails spectacularly in the chromosphere (ICarlsson & Steinl2002l) . 
Thus, these 3D simulations properly model the shock compres- 
sion and post-shock expansion, but inaccurately incorporate ra- 
diative heating, and give erroneous temperatures from the mass 
density and internal energy because of the assumption of Saha 
ionization equilibrium. 
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The ID simulations performed by ICarlsson & Steinl (Il997l) 
modeled the radiative losses in great detail, including the effect 
of slow hydrogen ionization/recombination on the equation of 
state. However, the ID geometry affected the temperatures in 
the chromosphere, leading to too hig h temperatures in the sh ocks 
as well as in the inter-shock phases. Leenaarts et al.l (1 20071) per- 
formed a 2D radiation-magneto-hydrodynamic (RMHD) simu- 
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lation with an EOS that took non-equilibrium hydrogen ioniza- 
tion into account using approximations for the radiation field 
based on I Solium! d 19991) . This simulation, however, still com- 
puted the radiative losses based on Saha equilibrium, and re- 
quired an ad-hoc heating term that prevented the temperature to 
drop below 2400 K. 

All of the above simulations were limited in such a way that 
an accurate prediction of the minimum temperature occurring in 
the quiet solar chromosphere could not be made. 

In this paper we discuss the temperature structure in a simu- 
lation that tries to remedy the various limitations of the previous 
models: it is two-dimensional, includes a non-equilibrium equa- 
tion of state, includes radiative losses without the assumption 
of Saha equilibrium and includes heating through absorption of 
coronal radiation. We do not include a significant magnetic field, 
so that our simulation serves as a baseline against which simula- 
tions with stronger field can be compared. 

In Sec. |2] we discuss the RMHD model assumptions in de- 
tail and discuss some of the basic physical processes that occur 
in the quiet chromosphere. In Sec. [3] we discuss the results of 
our run, especially paying attention to the accuracy of the radia- 
tive heating. We finish with a discussion and our conclusions in 
Sec.0 



2. The model 

We model the chrom osphere using the RMHD code Bifrost 
dGudiksen et al.ll2011l) . 



MHD equations. The Bifrost code solves the equations of resis- 
tive MHD including the effects of heat conduction and radiation: 



dp 

-£ + v • (pu) = o, 

ot 



(1) 

+ V • {puii - z) = -VP +jxB+pg, (2) 
+ V ■ (eu) = -P(V • a) - V • F c - V • F r + Q vhc + Q ]oule , (3) 



dpu 
~dt 

de 
di 

fi Q j = V x B, 
dB 

— = V x (a x B) - V x T/j, 
of = 



(4) 
(5) 



with p the mass density, u the velocity, x the viscous stress ten- 
sor, P the gas pressure, / the current density, B the magnetic 
field, g the gravitational acceleration, e the internal energy per 
volume, F c the energy flux due to heat conduction, F t the radia- 
tive energy flux, g v i sc the viscous energy dissipation, Qjouie the 
Joule heating, po the permeability of free space and 77 the electri- 
cal resistivity tensor. Bifrost employs a staggered grid and uses 
sixth-order operators to compute spatial derivatives. The time- 
ste pping is done using a third-order predictor-corrector s cheme 



stepping is done using a tnira-oraer preai 
bv lHvmanl d 1979b . modified for variable ti 



timestep. 



Equation of state. Conservation of internal energy is expressed 
as 



3kT 



e = 



n e + n H 2 + 



"1 

^ «,■ + «o 

1=1 



+ n m (e m +xm) 



"I 

+ Zj n iXi + e a (n e ,T). 



(6) 



Here, k, T, n e , nm, n\, n Q are Boltzmann's constant, the gas 
temperature, the electron density, the H2 molecule density, the 
number of levels in the hydrogen model atom, the density of 
atomic hydrogen in excitation or ionization state ; and the num- 
ber density of all other atoms and molecules that are not, or 
do not contain, hydrogen, respectively. The rotational and vi- 
brational energy per H2 molecule is em, Xm is the dissociation 
energy of H2, Xi I s me excitation or ionization energy of atomic 
hydrogen and e Q is the internal energy of all other atoms and 
molecules that are not, or do not contain, hydrogen. The last 
quantity depends on the temperature and electron density and 
is computed in LTE. 

Eq.|6]is solved together with the equations for charge conser- 
vation, hydrogen nucleus conservation, evolution equations for 
the atomic hydrogen level populations «, (5 bound levels plus 
the continuum): 



drij 
~dt 



V • (nm) = J] njP jt - m 2 Pij, (7) 

and an equation for non-equilibrium formation of H2 molecules: 
dn m 



dl 



+ V ■ («H2«) = Pf«i - Pb"l"H2, 



(8) 



where n\ is the ground state population of atomic hydrogen. 

The radiative part of the rat e coefficients P ,;, is computed us- 
ing the approximations given bv lSolluml (1 1 9991) . the temperature- 
dependen t rate coefficients Rf and R\, are taken from the UMIST 
database (IWoodall et al.ll2007l Iwww . udf a . netl l. 

These equations yield values for T, n e , nm, and n,. The gas 
pressure is then given by 



P = kT 



«e + «H2 + 



"I 



+ n. 



(9) 



and used in the momentum and energy equations (Eqs. |2H3}. 
This non-equilibrium equation-of-state requires advection of 6 
atomic and 1 molecular hydrogen population in addition to 
the 8 MHD variab l es, fo r a total of 15 adve cted quantities. 
See Leenaarts et al.l (120071) , iGoldind (1201 Ol) and lGudiksen et al.l 
(1201 lh for more details. 



Radiative heating. Formally, the radiative flux divergence is 
given by 



V-F, 



= I I a(v, n 
Jo Jn 



) (S(v,n) -I(v,n)) dfldv, 



(10) 



with 5, a and / the source function, opacity and intensity at 
frequency v in direction h into the solid angle D.. In practice, 
this double integral is too computationally expensive to com- 
pute, and various approximations are made. 

The integral over frequency is replaced by four representa- 
tive radiation bins, with their own associated bin-averaged opac- 
ity per mass unit /c,, photon destruction probability e, and thermal 
emission assumin g isotropic scattering and isotrop ic source 
function and opacity (lNordlundlll982t ISkartherJl2000l) . The ra- 
diative heating of such a so-called multi-group scheme is then 
given by 



Grad = 47Tp 2^ Ki(eiJi - E,) 



(11) 



;=i 



The implementati on of this scheme in the Bifrost code is dis- 
cussed in detail bv lHavek et al.l J2010h . The quantities a-,, e, and 
Ei are precomputed and tabulated assuming LTE populations 
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and line-scattering in the 2-level Van Regemorter approxima- 
tion. They are tabulated as function of the electron density and 
gas temperature . This is an improvement over previous work by 
iLeenaarts et al.l (120071) who used tables as function of the inter- 
nal energy and the mass density, implicitly assuming LTE values 
for the electron density and the temperature. 

The multi-group scheme has the drawback that to properly 
approximate the effects of the strongest spectral lines one needs 
to add several bins just to cover the largest opacities. Otherwise 
the cooling and heating effects of the few strong spectral lines in 
a bin otherwise comprised of spectral features with much lower 
opacity are washed out in the constructio n of the bin-average d 
radiation quantities (see e.g., Eq. 9-11 of lHavek et alJ (120101) '). 
Furthermore, the assumptions of an opacity in LTE and a source 
function from the 2-level Van Regemorter approximation do not 
work for the strongest chromospheric lines. In the case of simu- 
lations of the chromosphere this means that the radiative heating 
and cooling from the mid-chromosphere upwards is badly mod- 
eled with a multi-group scheme only. 

Therefore, the MHD model includes additional cooling in 
such lines in the following way: 

Optically thin radiative losses from the corona and transition 
region are included through a frequency-integrated loss function 
A(T) based on the coronal approximation for the level popula- 
tions: 

Gthin = -A(7> e n H . (12) 

The function A(T) is only significantly different from zero above 
T = 15000 K. 

Optically thick parameterized radiative losses and gains from 
H i, Mg n and Ca n are included through: 



G[H,Ca,Mg] - C(r)[ H ,Ca,Mg] «eP mhl 



111 



0.3 



(13) 



Here the constant k and the temperature-dependent coefficient 
C are determined from de tailed radiative transfer computations 
with t he RADYN (see e.gJCarlsson & Steinll995l) and Multi3D 
codes dLeenaarts & Carlssonll2009l) . and m. is the column mass. 
These functions include hydrogen lines and the Lyman contin- 
uum and the pertinent lines and continua from Ca n and Mg n 
(ICarlsson & Leenaartsll201 fl in preparation). 

Radiative heating through absorption of coronal radiation in 
UV continua is modeled through the representative bound-free 
absorption cross-section cr of He I at 25 nm: 

£>He = 0"He,25nm «He I e~ THc ' 25 " m /thin (14) 

with /thin the angle-averaged radiation field resulting from gthin 
dCarlsson & Leenaartsll201 fl in preparation). 

In addition we include an ad-hoc heating term that drives the 
temperature up to Tq = 1660 K once it drops below that value. 
It is given by: 

Q w = C w n 2 H {max(0,(T -T))] 2 , (15) 

with C w a constant that sets the heating rate. The simulation thus 
still has an artificially set mini mum temperature , but it is 740 K 
lower than in the simulation of ILeenaarts et al.l (120071) . We add 
this ad-hoc term to prevent our simulation from cooling to too 
low temperatures where the radiation tables we employ become 
inaccurate. 

All the different contributions are added so that the radiative 
heating term in Eq. [3]becomes 

- V • F t = Q mA + e, hin + Q H + Gca + QMg + fiHe + 2w (16) 
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Fig. 1. Snapshot of the simulated atmosphere. Top panel: gas 
temperature, with the column shown in Figs. [2] [3] and [7] — [T2l 
indicated by a dashed line. Bottom panel: -V ■ u, positive values 
indicate the gas is locally compressed in the co-moving frame. 

Simulation setup. We performed a 2D simulation with a grid 
size of 512 x 325, spanning from 1.5 Mm below (T500) = 1 to 
14 Mm above it in the lower corona. The horizontal grid spacing 
is 32.5 km, the vertical grid spacing is 28 km from the convection 
zone up to the low corona, and increased to 150 km higher up in 
the corona. 

The lower boundary is open, allowing fluid to freely enter 
and leave the box, while specifying the entropy of the inflow- 
ing gas to maintain sufficient energy flux into the computational 
domain. The upper boundary uses the methods of characteristics 
to extrapolate the MHD variables, letting waves exit the domain 
with almost no reflections. The upper boundary is set to strive 
to a temperature of 250,000 K to prevent the corona from cool- 
ing as a 2D weakly-magnetic simulation cannot sustain a corona 
self-consistently. We choose this rather low temperature because 
we model a weakly magnetic part of the atmosphere. Its exact 
value has no influence on the behaviour of the chromosphere. 
As formulated in the code, the thermal conduction operator re- 
quires a magnetic field to be present, so a weak magnetic field 
(average unsigned flux density in the photosphere of 0.3 G) is 
added, but is too weak to have any effect on the atmosphere. 

The simulation was started from a previously relaxed snap- 
shot computed without non-equilibrium hydrogen ionization and 
ran for 1 hour of solar time. We discard the first 10 minutes to 
remove start-up effects. 

3. Results 

Figure Q] shows the chromosphere in a snapshot from the sim- 
ulation after 34 min of solar time have elapsed. The upper 
panel shows the gas temperature, with granules visible below 
z = Mm, the high-temperature transition region and corona 
in white at the top. In between lies the chromosphere, visible 
as a cool background state pervaded by waves and shock-fronts 
with peak temperatures increasing with height. The bottom panel 
shows -V • it . a measure for the local compression rate of the 
gas. Intergranular lanes and their extension into the upper photo- 
sphere appear as weakly compressing purple stalks at the bottom 
of the panel. The scene is dominated by strongly compressive 
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Fig. 3. Energy balance along the cut indicated in Fig.Q] The up- 
per part shows heating, the lower part cooling. Black: total heat- 
ing rate (right-hand side of Eq. [3] excluding the ad hoc term gw- 
Green: compression work -P(V • u). Red: total radiative heating 
(Eq.fToli, excluding gw- Blue: ad-hoc heating gw (Eq.[T5Tl. 

shocks in the chromosphere and corona with expanding shock- 
wakes, such as the one around (x, z) = (13,2) Mm. Note that 
some of the shock fronts are co-spatial with the chromosphere- 
corona interface, and push the corona upward. An extreme ex- 
ample of this is seen along the dashed line at x = 13.1 Mm in 
the upper panel. The shock front at z = 3.7 Mm has pushed 
the corona up, leaving a large pocket of cool expanding gas in 
its wake. At z — 0.9 Mm a new shock front is propagating up 
through the cold gas, compressing it again. 

Figure |2] shows various quantities along the cut indicated by 
the dashed line in Fig. Q] Panel a shows the density profile. It 
decreases with height everywhere except around z = 3.8 Mm 
at the site of the shock front that pushes the corona upward. 
Panel b shows the total internal energy in black with its dis- 
tribution over various contributions. It has a peak at the shock 
front at z — 0.9 Mm and a smooth increase with height above 
it. The main contributors to the internal energy are the energy of 
the random motions of the gas particles (blue) and the ioniza- 
tion energy of hydrogen (red). The latter is the largest contribu- 
tion above 1 .6 Mm. There is a small but important contribution 
of H2 molecules at 0.5 Mm height. Panel c shows the vertical 
velocity. It exhibits a rough sawtooth shape common to shock- 
propagation in a stratified atmosphere. Panel d shows the tem- 
perature, with a high-temperature shock front at 0.9 Mm, and 
a plateau at 1660 K between 1 and 2.5 Mm. Panel e shows the 
gas pressure and panel f shows the compression rate, again indi- 
cating the presence of the two shock fronts at 0.9 and 3.8 Mm, 
weak compression at 2.8 Mm and expansion in the rest of the 
chromosphere. 

Figure [3] displays the energy balance in the chromosphere 
(the right-hand-side of Eq. [3} for the same column as Fig. [2] 
Viscous energy dissipation, Joule heating and thermal conduc- 
tion are all negligible in the chromosphere and are not shown. 
The energy balance is set by compression work and radiation. 




t [min] 

Fig. 4. Time evolution of the temperature at different heights in 
the column indicated by the dashed line in Fig. Q] The time mo- 
ment shown in that figure is indicated by the long tick marks at 
t - 8.4 min. Blue: z = 1 Mm; green: 1.5 Mm; red: 2 Mm. 
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Fig. 5. Diagram of the temperature occurrence as function of 
atmospheric height in the simulation. The red curve shows the 
temperature in the column marked in Fig. Q] The blue curves in- 
dicate the minimum and maximum temperatures as a function of 
height during the whole simulation run. The black curves shows 
the location of 50% He 1 (top curve) and H (bottom curve) ion- 
ization assuming Saha equilibrium and the average run of the 
electron density with height. 



As expected, the shock fronts are heated compressively and cool 
radiatively. The cool area between 1 and 3.4 Mm is cooling 
through expansion and is heated by radiation except at 2.7 Mm, 
where a horizontally propagating wave causes some compres- 
sion. The expansion cooling is stronger than the radiation heat- 
ing by an order of magnitude. This imbalance ultimately leads 
to the ad-hoc Qw contribution becoming active to prevent the at- 
mosphere from cooling below 1660 K. The instant in time shown 
in the figure shows a scene with large ad-hoc heating. However, 
this term is only intermittently active and switches off as soon as 
the temperature is above 1660 K. 

Figure |4] shows the time evolution of the temperature at dif- 
ferent heights along the dashed line of Fig. [TJ The blue curve 
between t — 16 and t = 24 minutes show a clean example of the 
temperature variations of passing shocks, with a rapid increase 
in temperature and a more gradual cooling phase after passage 
of the shock front. In general the time evolution is more irregu- 
lar because of slanted shock propagation and interference. The 
corona intermittently dips down to below 1.5 Mm, indicated by 
the red and green curves running off the temperature scale. 

The above discussion and Figs. |2]-0 again confirm the pic- 
ture of the quiet chromosphere as a layer that undergoes quick 



J. Leenaarts et al.: On the minimum temperature of the quiet solar chromosphere. 



5 





2 3 
z [Mm] 

Fig. 2. Properties of the atmosphere along the cut indicated in Fig. Q] a: Mass density; b: The solid black curve indicates the internal 
energy per mass unit, with the zero point at neutral atoms in the ground state. Red: ionization energy of hydrogen; black dotted: 
excitation energy of atomic hydrogen; blue: kinetic contribution of all atoms, molecules and electrons; green: contribution of the 
rotational, vibrational and dissociation energy of Ht molecules, c: Vertical gas velocity, positive means upward, d: Gas temperature, 
e: Gas pressure, f: compression rate -V • u. 



compression during shock passages followed by a longer phase 
of expansion cooling. The time scale of this expansion cooling 
tpdv = <?/(P(V • «)) varies from 200 s to 500 s assuming typical 
values of e = 10 12 erg g" 1 and 2 x 10 9 < P(V ■ u) < 5 x 10 9 
erg s _1 g -1 . Here e is roughly the thermal energy of solar gas 
at 12 000 K. The ionization energy remains nearly constant, so 
the gas in the chromosphere behaves approximately as an ideal 
gas. So, expansion cooling alone can, in cases of strong expan- 
sion, cool chromospheric gas to very low temperatures between 
2 shock passages, as the radiative heating rarely exceeds 2 x 10 8 
erg s _1 g _1 in the chromosphere. 

This rough estimate is confirmed by Fig. [5] It shows a his- 
togram of the occurrence of gas temperature values as a func- 
tion of atmospheric height. The range of temperatures increases 
with increasing height when going up from the photosphere. The 
maximum temperature curve shows the increase of peak shock 
temperature with height up to 1 Mm. The maximum temperature 
at 1 Mm is above 15 000 K, indicating that the corona occasion- 
ally reaches this far down. The thick dark band at 10000 K be- 
tween 1 and 4 Mm is caused by a combination of shock fronts 
and the layer of 10 000 K gas just below the corona (see the up- 
per panel of Fig.[TJ. This band is discussed in Sec. 13.21 Between 
z = 0.8 and 2.2 Mm the histogram peaks along a narrow dark 
band below T « 2000 K, indicating that such low temperatures 
are common in the simulated chromosphere. Above 2.2 Mm 
such low temperatures occur less frequently, but the atmosphere 
can be as cold as the ad-hoc heating threshold temperature of 
1660 K up to 3.4 Mm height. 

In the following subsections we investigate the accuracy of 
this result by discussing the various physical processes that de- 
termine the minimum temperature in the chromosphere and the 
accuracy with which they are represented in the simulation. 

3.1. H 2 thermostat action 

When the low-chromospheric temperature reaches down to 
2000-3000 K, H2 molecules can form in significant amounts. 
This process is exothermic, releasing 4.48 eV per molecule 
formed, compared to a thermal energy of kT — 0.19 eV per par- 
ticle at 2000 K. Thus, as H2 begins to form it releases a large 
amount of energy, which by Eq.|6]is predominantly converted to 
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Fig. 6. Diagram of the occurrence of temperature as function of 
mass density in the simulation. The red contours show the frac- 
tion of hydrogen atoms bound in H2 molecules, 2«h2/(«hi + 
2nm), as function of temperature and mass density, assuming 
ICE and all hydrogen neutral. The contours are spaced a factor 
100 apart, with the labels indicating the exponent a in 10 fl . The 
blue line at T = 1 .66 kK specifies the threshold for the ad hoc 
heating (Eq.[T5l). 



thermal energy, counteracting the expansion cooling. This effect 
is illustrated in Fig. [6] It shows the simultaneous occurrence rate 
of combinations of temperature and mass density in the simula- 
tion, with overplotted contours of the fraction of hydrogen atoms 
bound in H2 assuming instantaneous chemical equilibrium. 

The minimum temperature in the chromosphere follows the 
curve that represents 10% of all hydrogen bound in H2 between 
p = 10~ 9 g cm~ 3 and 10~ 75 g cm 4 . At lower mass densities 
(higher up in the chromosphere) the formation rate of H2 be- 
comes so low that it cannot form in large enough amount to pre- 
vent the expanding parts of the atmosphere from cooling further. 
This is indicated by the turnoff of the bottom of the grey cloud 
away from the red curve at p = 10~ 9 g cirT 3 . Once the temper- 
ature drops below 1660 K the ad-hoc heating becomes active, 
preventing the atmosphere from cooling even further. 
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Fig. 7. Various contributions to the total radiative heating rate 
(Eq.fToli for the column indicated in Fig.Q] Black solid: total ra- 
diative heating. Black dashed: Radiative heating from the multi- 
group scheme (Q m d, see Eq. ITTb . Red: Qh- Green: <2ca- Blue: 
2thin- Red dashed: g He - Grey with black dots: Q Ms . 



3.2. H and He as thermostats 

Figure Q] shows a layer of 10000 K hanging as a ragged skirt 
below the corona. This layer is one of the causes of the dark 
band at the same temperature in Fig. The overlaid He ion- 
ization curve demonstrates that Hei works similarly to H2 as 
a thermostat in this simulation. This is unphysical, however, 
since the actual ionization-recombination balancing is likely to 
be even more out of instantaneous equilibrium for helium than 
in the case of hydrogen. Thus, the Saha equilibrium assumed in 
our simulation for helium is erroneous. Relaxing this assump- 
tion is likely to remove helium's ther mostat action, just as for 
non-equilibrium hyd rogen ionization (ICarlsson & Steinl l2002t 
Leenaarts et al .120071) . Indeed, hydro gen does not cause any such 
thermostat clustering of temperature values (which would be at 
the lower black curve in Fig. [5} since it is treated properly in 
non-equilibrium. The erroneous assumption of LTE balancing 
for He 1 does not affect the present analysis, however, since it af- 
fects only the hottest phases of the chromospheric gas and not 
the coolest ones addressed here. 

3.3. Radiative heating 

The chromosphere is heated by radiation during its cool phases. 
Fig.Qshows the various contributions that make up the total ra- 
diative heating rate in our simulation along the column indicated 
in Fig.Q] 

The radiative heating due to the multi-group scheme (Eq.QT} 
dominates the heating below z — 1 Mm as indicated by the 
near-equality of the solid and dashed black curves. It also cools 
strongly in the upper chromosphere. Absorption of coronal ra- 
diation (red dashed curve, Eq. fl4t adds heating above 1.3 Mm 
and is dominant between 1.8 and 3.2 Mm. The Can lines (green, 
Eq. [T3l cool in the shock front at 0.9 Mm and above 2.8 Mm, 
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Fig. 8. Comparison of the multi-group radiative heating rate in 
the first bin (black solid, Q\ see Eq.fTTI) compared to the heating 
rate due to H in NLTE (red) and LTE (green). 



they are the dominant contribution in the cool gas between 1 
and 1.8 Mm. Hydrogen (red, Eq.QjJi cools in the shock front at 
0.9 Mm and above 2.8 Mm, where the Lyman lines and continua 
become effectively optically thin. Optically thin losses from the 
corona (blue) do not play a role in the chromosphere. The Mg 11 
(grey-black dotted) lines only cool significantly in the upper 
chromosphere. 

We will now discuss the different contributions in detail and 
show them for the representative column from our simulation 
displayed in Figs.|2]-|4] 

Negative hydrogen ion H . Absorption of radiation in the H 
continuum is a potentially large source of heating. This absorp- 
tion is included in the first radiation bin of our model. We com- 
puted the H heating in detail for the snapshot of Fig.[T]assuming 
LTE, i.e., the source function is the Planck function and the H 
density is set by its Saha-equilibrium with neutral hydrogen. For 
comparison, we computed the H heating rate in NLTE includ- 
ing the following reactions: 

H + hv <-> H + e, 
H+H <-> 2H + e, 

H + e <+> H + 2e, 

H+p^ 2H, 
H+H <-> H 2 + e, 

where we kept the neutral hydrogen, proton, electron and, 
H2 density constant. This constancy is justified because of 
the high number density of these particles relative to the 
H d ensity. The reaction rate s for the latter four reactions are 
from lLambert & Pagell (119681) . 

Figure [8] shows a comparison with the first bin of the multi- 
group heating rate (black, Eq.QT]and the heating rate caused by 
H bound-free and free-free transitions, assuming LTE (green) 
and NLTE (red). The LTE assumption generates the largest heat- 
ing, caused by both the low Planck function compared to the av- 
erage radiation field, and the high H number density. The NLTE 
heating rate is much smaller. This is caused by a combination 
of photo-ionization of H , leading to a NLTE departure coeffi- 
cient much smaller than unity, and strong scattering, leading to 
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Fig. 9. Comparison of the parameterized radiative heating rate 
due to Can (Qcn, see Eq. [L31 to the heating computed based 
on a detailed radiative transfer computation. Black: Qc d - Blue: 
heating due to Can lines. Green: heating due to Can continua. 
Red: total Can heating. The blue curve is nearly equal to the red 
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Fig. 10. Comparison of the parameterized radiative heating rate 
due to Mgn (2ivig, see Eq. [T3l ) to the losses computed based 
on a detailed radiative transfer computation. Black: Qwig- Blue: 
heating due to Mg n lines. Green: heating due to Mg n continua. 
Red: total Mg n heating. The blue curve is nearly equal to the red 



smaller J v -S v splitting. Absorption of radiation by H is there- 
fore not efficient in heating the chromosphere. 

The first bin in the multi-group scheme closely matches 
the heating due to H where H dominates the opacity (up to 
0.8 Mm). It deviates in the shock around 0.9 Mm and above. 
The first bin underestimates the NLTE H heating in the cool 
area between 1 and 2.5 Mm. Above 3 Mm it erroneously over- 
estimates the NLTE cooling by several orders of magnitude. 

Ca n. The lines of Ca n are other strong heating agents in the 
chromosphere. We computed the Ca n heating rates for the snap- 
shot of Fig. Q] using the radiative transfer code Multi3d. This 
detailed computation was performed treating each column in the 
simulation as a plane-parallel atmosphere and included the ef- 
fects of partial redistribution (PRD) in the H& K lines. The re- 
sults for the representative column are shown in Fig. |9] The Ca n 
continua (green curve) have a negligible effect on the total heat- 
ing rate. The lines cool strongly in the shock front at 0.8 Mm 
and above 3.4 Mm. They heat the atmosphere in the cool area 
between 0.9 and 3.4 Mm. 

The parameterized Ca n cooling employed in our simulation 
is in reasonable agreement with the detailed computation above 
0.5 Mm height. The largest differences are the overestimation of 
the cooling between 3.5 and 4.2 Mm and shifted location around 
3 Mm where the heating term switches sign. The difference be- 
low 0.5 Mm is caused by a cutoff to avoid doubling the heating 
rate as heating at low heights is also included in the multi-group 
scheme. 

Manual inspection of the heating rate in many different 
columns of the simulation shows that the parameterized heating 
is accurate most of the time. 

Mg ii. Another potentially large contributor to the total chromo- 
spheric heating is Mg n. We computed the radiative losses in de- 
tail including the effects of PRD in the h&k lines. The result 
for the representative column is given in Fig. [10] The h&k lines 
cool in the shock front at 0.9 Mm, heat the cool area above the 
shock front and cool the upper chromosphere. The continua play 
a minor role. 
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Fig. 11. Comparison of the parameterized radiative heating rate 
due to hydrogen (Qh, see Eq. [T3l to the losses implied by the 
non-equilibrium hydrogen radiative rates. Black: Qu- Blue: im- 
plied heating due to hydrogen lines. Green: implied heating due 
to hydrogen continua. Red: total implied heating. 



The cutoff in the parameterized heating pushes the heating 
to zero at 1.5 Mm. The heating is much lower than the de- 
tailed computation in the cool area between 1 .5 and 3 Mm and it 
shows the same shift as for Ca n in the location where the heating 
switches sign. The cooling peak around 3.8 Mm is reproduced 
correctly by Q M „. 

Manual inspection shows that the parameterized Mg n heat- 
ing is accurate for gas temperatures above 6000 K. Below this 
temperature the heating rate is modeled qualitatively correct, but 
can deviate because the effects of the precise atmospheric struc- 
ture, in particular the velocity field are not taken into account in 
the simple parameterization. 

Hydrogen. Hydrogen has only a small effect on the heating rate 
of the lower and middle chromosphere in the semi-empirical 
time- independent VAL C model atmosphere dVernazza et al.l 
1 198 lh . In this model the Lya transition is so optically thick that 
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it does not play a role in the energy balance, and the losses in Ha 
are approximately offset by the absorption in the Balmer contin- 
uum. This effect remain s valid in the time-dependent models of 
ICarlsson & Steinl d2002l) computed with the RADYN code, ex- 
cept in shock fronts. Unfortunately it is currently not possible to 
compute the detailed time-dependent heating rate due to hydro- 
gen in multi-dimensional simulations, so we cannot rigorously 
prove the same effect holds in our 2D simulation. However, the 
small hydrogen heating rate is mostly an effect of the atomic 
structure of hydrogen, and does not depend strongly on the exact 
chromospheric structure. Therefore, our 2D simulation should 
behave in a similar manner. This is confirmed in Fig. QT] 

This figure shows the parameterized heating due to hydrogen 
(Eq. IT3l> in black. The shock front at 0.9 Mm is cooling, there is 
little heating between 1.2 and 2.5 Mm, and slowly increasing 
cooling above 2.8 Mm, peaking at 4.1 Mm where the Lya line 
becomes optically thin. 

Figure. QT| also shows the heating implied by the radiative 
transition rates computed from Eq. [7] We define the implied 
heating rate in a transition between levels ; and j (i < j) by 

Gimp =hv (Rij-Rji), (17) 

with vo the line center or ionization edge frequency, and R,j 
the radiative rate from level i to level j. This simple defini- 
tion makes an error for bound-free transitions due to the non- 
negligible width of the ionization edges, but this error is small. 
Note that the implied heating does not appear as a source term 
in Eq. [161 the assumed radiation field only serves to define the 
radiative rates in Eq. [7] It reproduces the near-zero heating in the 
cool area between 1 Mm and 2.8 Mm. The shock-front and the 
upper chromosphere are heated. 

We conclude that both the parameterized cooling and the im- 
plied heating rate reproduce the lack of hydrogen heating in the 
cool phases of the chromosphere. Any errors caused be the ap- 
proximations have negligible effect on the minimum temperature 
in our model due to the small heating rate compared to other 
heating agents. 

Absorption of coronal radiation. Nearly all coronal radiation 
emitted towards the chromosphere is eventually absorbed (a 
small part is backscattered into space). Exactly where this ra- 
diation is absorbed depends on its detailed spectral energy dis- 
tribution and the absorption coefficient in the chromosphere. 
However, these details are of minor importance to the overall 
heating rate as long as all energy is absorbed. Our simulation 
correctly models this behavior. The finite extent of the corona in 
the simulation and its relatively low temperature might result in a 
too low amount of radiation absorbed in the chromosphere. The 
coronal loss function {Q^ m , Eq.[T2l) peaks strongly in the lower 
corona due to its quadratic density dependence, so we expect this 
error to be small. 

4. Discussion and conclusions 

Robustness of the minimum temperature The formation of H 2 
acts as a thermostat in the lower chromosphere, preventing the 
temperature from dropping below 2000 K. ICE is not valid in the 
middle and upper chromosphere. There, the slow formation rate 
of H2 prevents thermostat action, allowing solar gas to cool well 
below 2000 K. Our simulation correctly models this behavior. 

In Sec. 13.31 we discussed the various lines and bound- 
free transitions that can radiatively heat the chromosphere. We 
showed heating in lines of Ca 11 and absorption of coronal radi- 
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Fig. 12. Molecular densities in the exemplary column indicated 
by a dashed line in Fig.Q] The densities of CO (green), OH (blue) 
and TiO (red) are based on instantaneous chemical equilibrium. 
The H2 density (black dashed) is based on our EOS. The total 
number density of hydrogen atoms is plotted in solid black for 
comparison. 



ation are the dominant processes. We showed that H is inef- 
ficient in heating cool pockets of chromospheric gas owing to 
strong scattering and low H population due to photo-ionization. 
Mg 11 and hydrogen play a minor role in the heating of the chro- 
mosphere. Our simulation accurately models the heating due to 
Can and coronal radiation. The other processes are less accu- 
rately modeled but these have, in sum, only a small effect on the 
energy balance. 

The heating due to acoustic waves is self-consistently in- 
cluded in the simulation since we include the upper convection 
zone and the stochastic excitation of acoustic waves there. The 
limited resolution may lead to an underestimate of the excitation 
of high-frequency waves. A wave at a frequency of 20 mHz is 
represented with 1 1 grid-points so the effect will only be impor- 
tant for really high frequency waves where simulations and ob- 
servations i ndicate the effect for the heating of the chromosp here 
is minimal dFossum & Carlssonll2005t [Carlsson et al.ll2.Q07D . 

The ad-hoc heating term acts intermittently in the largest 
expansion bubbles, thus keeping our minimum chromospheric 
temperature at 1660 K. 

We do not expect that a 3D simulation with similar physics 
will change this resul t. Three-dimens i onal s imulations, such as 
the one reported on bv lLeenaarts et al.l(l2010l) show temperatures 
down to 2000 K even with LTE hydrogen ionization and instan- 
taneous H2 formation. It is very likely that such 3D simulations 
with a non-equilibrium EOS would s how the same low tem per- 
atures as reported here (see Fig. 4 of Leenaarts et al.l (120071) for 
an illustration of the resulting temperature differences between 
these equations-of-state). 

We therefore conclude that our simulation provides an upper 
bound to the minimum temperature of the non-magnetic chro- 
mosphere. The simulation predicts the occurrence of gas tem- 
peratures as low as 1660 K up to 3.4 Mm height. This low tem- 
perature is common between 0.7 and 2.5 Mm height. Our ad-hoc 
heating prevents our model from cooling further. This means that 
a non-magnetic chromosphere can have temperatures even lower 
than this value. It seems impossible to avoid such low tempera- 
tures using only radiation and hydrodynamic processes. 
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Is the chromosphere really this cool? The sun has a tur- 
bulent photospheric magn etic field of the order of 130 G 
dTruiillo Bueno et aUl2004l) . This turbulent field is likely to ex- 
tend into the quiet chromosphere and might generate enough 
Joule heating to prevent the quiet chromosphere from cooling 
down to temperatures as low as in our simulation. In addition the 
presence of a relatively strong magnetic field will also change 
the propagation properties of the generated acoustic waves, re- 
ducing their amplitudes by confining them to flux tubes and/or 
by conversion to other magneto-hydrodynamic modes in the 
vicinity of the /3 = 1 surface. 

The simulation s of a s olar surface dyna m o by 
IVogler & Schiisslerf (120071) and ISchussler & Voderl (120081) 
suggest that the strength of the turbulent field declines rapidly 
with height, which would limit the amount of mid and high 
chromosphe ric heating the fiel d could cause. On the other hand, 
the work of ISchriiver & Tiflel (120031) suggests that part of the 
network magnetic fields can connect back in the internetwork 
photosphere. Such a multiscale magnetic carpet would lead 
to an increase of the field strength in the quiet chromosphere 
relative to a purely locally generated field, with a corresponding 
increase in the potential for Joule heating. 

Magnetic fields in the low ionization-degree chromosphere 
may also lead t o the generation of electri c currents through 
neutral-ion drag (iKrasnoselskikh et al.ll2010f) . 

Our simulation serves as a baseline case to test the effect of 
such magnetic heating proc esses. It should be compared against 
3D models like the one by IVogler & Schiissleil (120071) but ex- 
tended up into the corona, and including all the physics dis- 
cussed in this paper to study the effect of the turbulent field. 
Such a simulation does not require a large horizontal extent of 
the computational domain and, while computationally very de- 
manding, can in principle be done on the largest currently avail- 
able supercomputers. 

Simulation of the effect of the magnetic carpet requires a 
larger horizontal extent to include some network magnetic field 
and a large enough area of quiet sun. The high spatial resolution 

10 km grid spacing) needed to get local dynamo action com- 
bined with the large size of a network cell (« 30 Mm) makes 
such a simulation challenging. 

At this moment the minimum temperature of the quiet sun 
chromosphere remains an open question. Our simulation shows 
a non-magnetic chromosphere will get colder than our ad-hoc 
limit of 1660 K, but it is unknown whether regions uninfluenced 
by magnetic fields occur in the chromosphere. However, if these 
low temperature areas exist, they can in principle be observed. 

Suggestions for observations Direct observation of cool pock- 
ets of chromospheric gas is hard. Atomic spectral lines with 
sufficient opacity form generally in NLTE, making it hard 
to derive temperatures. The Atacama Large Millimeter Array 
(ALMA) would be the ideal instrument to probe chromo- 
spheric temperatures because of its high resolution, the LTE 
source function and relatively well-d efined formation height 
range of the sub-millimeter c ontinua (ILoukitcheva et al.ll2.Q04t 
IWedemever-Bohm et al.ll2007l) . 

Alternatively, off-limb spectroscopy in molecular lines with 
sub-arc-second resolution and low stray light could be used. 
As an example we show the number densities of some abun- 
dant molecules as a function of height in our exemplary col- 
umn in Fig. [12] Because the chemical reaction timescales be- 
come very large in the chromosphere the number densities 
for the species computed assuming ICE should be taken as a 



maximum possible value only, and are likely orders of mag- 
nitude too high. The molecules of choice for such observa- 
tions seem to be H 2 and CO. The sol ar spectrum shows H2 
UV emission lines dJordan et al.l 1 19781) that can be observed 
by the SUMER instrument aboar d the Solar and Heliosp heric 
Observatory (SOHO) spacecraft dKuhn & Morgan! 120061) . The 
mere presence of molecular lines at several arc-seconds above 
the limb would prove the existence of low temperature gas at 
high-chromospheric heights. Unfortunately the reverse is not 
true, the absence of lines might merely indicate the absence of 
molecules and not the absence of cool gas. 

Conclusions We have performed a 2D radiation-MHD simu- 
lation of the non-magnetic solar chromosphere, including the 
effects of non-equilibrium ionization of hydrogen and non- 
equilibrium formation of H2 molecules and NLTE radiative 
transfer. We find that our simulation contains pockets of cool gas 
with temperatures down to 1660 K from 0.8 Mm up to 3.4 Mm 
height. We discuss the physical mechanisms that set the mini- 
mum temperature and find that our simulation likely overesti- 
mates the minimum temperature. We conclude it is impossible 
to avoid such low temperatures with radiation-hydrodynamical 
processes only. 

Such cool pockets of chromospheric gas might be even 
cooler in the real sun than in our simulation, provided a quiet 
chromosphere without significant magnetic heating exists. We 
suggest off-limb molecular spectroscopy to look for such pock- 
ets, and suggest 3D simulations with sufficient resolution to sup- 
port a local dynamo with and without network-like magnetic 
fields to investigate whether the assumption of negligible mag- 
netic heating is justified. 
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